#include "transform.h"

#include <cmath>

#ifndef M_PI
#define M_PI		3.14159265358979323846
#endif

namespace raster {

using mathengine::Vec2f;
using mathengine::Vec3f;
using mathengine::Vec4f;
using mathengine::Matrix44f;

Transform operator*(float s, const Transform& t)
{
    return t * s;
}

/// T(tx,ty,tz) : P(x,y,z) => P'(x',y',z')
/// (x',y',z',1)^T = T * (x,y,z,1)^T
/// T =
/// (1 0 0 tx
///  0 1 0 ty
///  0 0 1 tz
///  0 0 0 1)
/// T^{-1} = T(-tx,-ty,-tz) =
/// (1 0 0 -tx
///  0 1 0 -ty
///  0 0 1 -tz
///  0 0 0 1)
/// \see https://www.cnblogs.com/fortunely/p/17749241.html
Transform Translate(float tx, float ty, float tz)
{
    Matrix44f m(
        {1, 0, 0, tx,
         0, 1, 0, ty,
         0, 0, 1, tz,
         0, 0, 0, 1});
    Matrix44f im(
        {1, 0, 0, -tx,
         0, 1, 0, -ty,
         0, 0, 1, -tz,
         0, 0, 0, 1});
    return Transform(m, im);
}

Transform Translate(const Vec3f &tv)
{
    return Translate(tv.x(), tv.y(), tv.z());
}

/// (x',y',z',1)^T = S * (x,y,z,1)^T
/// S =
/// (sx 0   0  0
///  0  sy  0  0
///  0  0  sz  0
///  0  0   0  1)
///
/// \see https://www.cnblogs.com/fortunely/p/17749241.html
Transform Scale(float sx, float sy, float sz)
{
    if (0 == sx || 0 == sy || 0 == sz) {
        return {};
    }
    Matrix44f m(
        {sx, 0,  0,  0,
         0,  sy, 0,  0,
         0,  0,  sz, 0,
         0,  0,  0,  1});
    Matrix44f im(
        {1/sx, 0,    0,  0,
         0,    1/sy, 0,  0,
         0,    0,    sz, 0,
         0,    0,    0,  1});
    return Transform(m, im);
}

Transform Scale(const Vec3f &sv)
{
    return Scale(sv.data[0], sv.data[1], sv.data[2]);
}

/// (x',y',z',1)^T = R * (x,y,z,1)^T
/// R(θ) =
/// (1 0    0     0
///  0 cosθ -sinθ 0
///  0 sinθ cosθ  0
///  0 0    0     1)
///
Transform RotateX(float theta)
{
    float rad = theta * M_PI / 180;
    Matrix44f m(
        {1, 0,               0,              0,
         0, std::cos(rad),   -std::sin(rad), 0,
         0, std::sin(rad),   std::cos(rad),  0,
         0, 0,               0,              1});
    Matrix44f im(
        {1, 0,               0,              0,
         0, std::cos(rad),   std::sin(rad),  0,
         0, -std::sin(rad),  std::cos(rad),  0,
         0, 0,               0,              1});
    return Transform(m, im);
}

/// (x',y',z',1)^T = R * (x,y,z,1)^T
/// R(θ) =
/// (cosθ  0  sinθ 0
///  0     1  0    0
///  -sinθ 0  cosθ 0
///  0     0  0    1)
Transform RotateY(float theta)
{
    float rad = theta * M_PI / 180;
    Matrix44f m(
        {std::cos(rad),  0,  std::sin(rad),  0,
         0,              1,  0,              0,
         -std::sin(rad), 0,  std::cos(rad),  0,
         0,              0,  0,              1});
    Matrix44f im(
        {std::cos(rad),  0,  -std::sin(rad), 0,
         0,              1,  0,              0,
         std::sin(rad),  0,  std::cos(rad),  0,
         0,              0,  0,              1});
    return Transform(m, im);
}

/// (x',y',z',1)^T = R * (x,y,z,1)^T
/// R(θ) =
/// (cosθ  -sinθ 0  0
///  sinθ  cosθ  0  0
///  0     0     1  0
///  0     0     0  1)
Transform RotateZ(float theta)
{
    float rad = theta * M_PI / 180;
    Matrix44f m(
        {std::cos(rad),  -std::sin(rad),  0,  0,
         std::sin(rad),  std::cos(rad),   0,  0,
         0,                0,             1,  0,
         0,                0,             0,  1});
    Matrix44f im(
        {std::cos(rad),  std::sin(rad),   0,  0,
         -std::sin(rad), std::cos(rad),   0,  0,
         0,                0,             1,  0,
         0,                0,             0,  1});
    return Transform(m, im);
}

Transform Rotate(float nx, float ny, float nz, float theta)
{
    return Rotate(Vec3f{nx, ny, nz}, theta);
}

/// (x',y',z',1)^T = R * (x,y,z,1)^T, rotate axis = n(nx, ny, nz), rotate angle = theta
/// R(n; theta) =
/// (nxnx(1-cosθ)+cosθ   nxny(1-cosθ)-nz*sinθ  nxnz(1-cosθ)+nysinθ   0
///  nxny(1-cosθ)+nzsinθ nyny(1-cosθ)+cosθ     nynz(1-cosθ)-nxsinθ   0
///  nxnz(1-cosθ)-nysinθ nynz(1-cosθ)+nxsinθ   nz(1-cosθ)+cosθ       0
///  0                   0                     0                     1)
///
/// \see https://www.cnblogs.com/fortunely/p/17749241.html
Transform Rotate(const Vec3f &axis, float theta)
{
    Vec3f n = axis.normalized();
    float nx = n.x(), ny = n.y(), nz = n.z();
    float rad = theta * M_PI / 180;
    float cos = std::cos(rad);
    float sin = std::sin(rad);

    Matrix44f m(
        {nx*nx*(1-cos)+cos,     nx*ny*(1-cos)-nz*sin,  nx*nz*(1-cos)+ny*sin,  0,
         nx*ny*(1-cos)+nz*sin,  ny*ny*(1-cos)+cos,     ny*nz*(1-cos)-nx*sin,  0,
         nx*nz*(1-cos)-ny*sin,  ny*nz*(1-cos)+nx*sin,  nz*(1-cos)+cos,        0,
         0,                     0,                     0,                     1});
    Matrix44f im(m.inverse());
    return Transform(m, im);
}

/// M = T * R * S
///
/// \see https://www.cnblogs.com/fortunely/p/18709389
Transform Modeling(float tx, float ty, float tz, float sx, float sy, float sz, float nx, float ny, float nz, float theta)
{
    return Translate(tx, ty, tz) * Rotate(nx, ny, nz, theta) * Scale(sx, sy, sz);
    // return Scale(sx, sy, sz) * Rotate(nx, ny, nz, theta) * Translate(tx, ty, tz); // for test @Martin
}

Transform Modeling(const Vec3f &translate, const Vec3f &scale, const Vec3f &axis, float theta)
{
    return Modeling(
        translate.x(),  translate.y(),  translate.z(),
        scale.x(),      scale.y(),      scale.z(),
        axis.x(),       axis.y(),       axis.z(),       theta);
}

/// Convert P(x,y,z) from world coordinate system to view coordinate system.
///
/// \see https://www.cnblogs.com/fortunely/p/18709389
Transform Camera(const Vec3f &eye_pos, const Vec3f &gaze_dir, const Vec3f &view_up)
{
    Vec3f w = -gaze_dir;
    Vec3f u = crossProduct(view_up, w);
    Vec3f v = crossProduct(w, u); // view_up

    return CoordAToCoordB(u, v, w, eye_pos);
}

/// P' = M_{ortho,norm}P
/// M_{ortho,norm} =
/// (2/(r-l) 0       0       -(r+l)/(r-l)
///  0       2/(t-b) 0       -(t+b)/(t-b)
///  0       0       2/(n-f) (n+f)/(n-f)
///  0       0       0       1)
///
/// \see https://www.cnblogs.com/fortunely/p/18709389
Transform Orthographic(float left, float right, float bottom, float top, float near, float far)
{
    float l = left, r = right, b = bottom, t = top, n = near, f = far;

    if (l == f || t == b || n == f) { return {}; }

    Matrix44f m(
        {2/(r-l),  0,        0,        -(r+l)/(r-l),
         0,        2/(t-b),  0,        -(t+b)/(t-b),
         0,        0,        2/(n-f),  (n+f)/(n-f),
         0,        0,        0,        1});
    return Transform(m, m.inverse());
}

Transform Orthographic(const Vec3f &lbn, const Vec3f &rtf)
{
    float left = lbn.x(), bottom = lbn.y(), near = lbn.z();
    float right = rtf.x(), top = rtf.y(), far = rtf.z();
    return Orthographic(left, right, bottom, top, near, far);
}

/// P' = M_{pers, norm}P
/// M_{pers, norm} =
/// (-2n/(r-l)  0          (r+l)/(r-l)   0
///  0          -2n/(t-b)  (t+b)/(t-b)   0
///  0          0          (n+f)/(n-f)   -2nf/(n-f)
///  0          0          -1            0)
///
/// \see https://www.cnblogs.com/fortunely/p/18709389
Transform Perspective(float left, float right, float bottom, float top, float near, float far)
{
    float l = left, r = right, b = bottom, t = top, n = near, f = far;
    if (l == f || t == b || n == f) { return {}; }

    Matrix44f m(
        {-2*n/(r-l),  0,           (r+l)/(r-l),  0,
         0,           -2*n/(t-b),  (t+b)/(t-b),  0,
         0,           0,           (n+f)/(n-f),  -2*n*f/(n-f),
         0,           0,           -1,           0});
    return Transform(m, m.inverse());
}

/// \brief Perspective transformation.
/// Assuming left + right = bottom + top = 0
///
/// \param near the distance between the near plane and the view origin
/// \param far  the distance between the far plane and the view origin
/// \param fovY Field of View Y, unit: degree
/// \param aspect the aspect ratio of the clipping window
///
/// \see https://www.cnblogs.com/fortunely/p/18709389
/// \see https://www.cnblogs.com/fortunely/p/18114286
Transform Perspective(float neardis, float fardis, float fovY, float aspect)
{
#if 0
    float n = -neardis, f = -fardis;
    float fovYRad = fovY * M_PI / 180;
    float t = std::tan(fovYRad / 2) * neardis;
    float b = -t;
    float r = aspect * t;
    float l = -r;

    return Perspective(l, r, b, t, n, f);
#else
    float n = -neardis, f = -fardis;
    float fovYRad = fovY * M_PI / 180;
    float cot = 1 / std::tan(fovYRad / 2);
    Matrix44f m(
        {cot/aspect, 0,    0,            0,
         0,          cot,  0,            0,
         0,          0,    (n+f)/(n-f),  -2*n*f/(n - f),
         0,          0,    -1,           0});
    return Transform(m, m.inverse());
#endif
}

/// \brief Coordinate system transform. From P(x,y,z)^A in coordinate system A[O(0,0,0);e1,e2,e3]
/// to P'(x',y',z')^B in coordinate system B[O'(x0,y0,z0);u,v,n].
///
/// M_{xyz,x'y'z'} =
/// (ux uy uz -u*O'
///  vx vy vz -v*O'
///  nx ny nz -n*O'
///  0  0  0  1)
///
/// \see https://www.cnblogs.com/fortunely/p/17761095.html
Transform CoordAToCoordB(const Vec3f &u, const Vec3f &v, const Vec3f &n, const Vec3f &b_origin)
{
    Vec3f uu = u.normalized();
    Vec3f vv = v.normalized();
    Vec3f nn = n.normalized();

    Matrix44f m(
        {uu.x(), uu.y(), uu.z(), -(uu * b_origin),
         vv.x(), vv.y(), vv.z(), -(vv * b_origin),
         nn.x(), nn.y(), nn.z(), -(nn * b_origin),
         0,      0,      0,      1});
    return Transform(m, m.inverse());
}

Transform Viewport(int width, int height)
{
    return Viewport(width, height, mathengine::Vec2f({0, 0}));
}

/// Transform object from CVV(canonical view volume) to viewport of the screen, considering pixel size.
/// Def. ignore pixel size
///
/// M_{vp} =
/// (nx/2 0     0  x0+nx/2
///  0    -ny/2  0  y0+ny/2
///  0    0     1  0
///  0    0     0  1)
///
/// \param width the width of the viewport rect, unit: pixel
/// \param height the height of the viewport rect, unit: pixel
/// \param left_corner the screen coordinate of the left corner of the viewport rect
///
/// \see https://www.cnblogs.com/fortunely/p/18106741
Transform Viewport(int width, int height, const mathengine::Vec2f& left_corner)
{
    float nx = width;
    float ny = height;
    float x0 = left_corner.x();
    float y0 = left_corner.y();
    bool ignorePixelSize = true;

    Matrix44f m;
    if (ignorePixelSize) {
        m = Matrix44f(
            {nx/2, 0,     0,  x0 + nx/2,
             0,    ny/-2,  0,  y0 + ny/2,
             0,    0,     1,  0,
             0,    0,     0,  1});
    } else {
        m = Matrix44f(
            {nx/2, 0,     0,  x0 + (nx-1)/2,
             0,    ny/-2,  0,  y0 + (ny-1)/2,
             0,    0,     1,  0,
             0,    0,     0,  1});
    }

    return Transform(m, m.inverse());
}

/// Vertex transform to normal transform of the vertex of the surface
///
/// \see https://www.cnblogs.com/fortunely/p/18178861
Transform ToNormalTransform(const Transform &vertexTransform)
{
    return vertexTransform.toNormal();
}

}
